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Summary of Progress Prior to This Period 

The main purpose of this research has been to develop a rigorous theory and corresponding 
computational algorithms for through-the-thickness analysis of composite plates. This type of 
analysis is needed in order to find the elastic stiffness constants for a plate and to post-process the 
resulting plate solution in order to find approximate three-dimensional displacement, strain, and 
stress distributions throughout the plate. 

We have settled on the variational-asymptotical method (VAM) 1 as a suitable framework in 
which to solve these types of problems. The VAM was applied to laminated plates with constant 
thickness in the work of Atilgan and Hodges 2 . The corresponding geometrically nonlinear global 
deformation analysis of plates was developed by Hodges, Atilgan, and Danielson 3 . A different 
application of VAM, along with numerical results, was obtained by Hodges, Lee, and Atilgan 4 . 
(Copies of these papers have been delivered to Mr. Hinnant.) 

In Ref. 2, the “first approximation” is exactly the same as classical laminated plate theory. 
The “second approximation” takes transverse shear deformation into account and is developed only 
for plates with certain restrictions in their construction. To remove the restrictions one must “kill” 
certain interaction terms in the strain energy, and the means for doing so for general laminated 
plates were not given in that paper. 

In Ref. 3, a set of kinematical and intrinsic equilibrium equations are derived for large deflec- 
tion and rotation but with small strain. The relationship between the drilling rotation and the 
other kinematical variables gives new insight into the drilling moment and its role in beam-plate 
connectivity. This work has shown that drilling type rotation is not an independent degree of 
freedom in plate theory. An applied drilling moment at a point on a plate is not resisted at all by 
the plate. Such a moment, in order to have any physical resistance from the plate, must be applied 
over a finite area. Other than this, a point drilling moment can only be resisted by a plate if the 
plate model is derived from couple-stress elasticity 5 . 

The development in Ref. 4 includes transverse shear in the “first approximation” and is stopped 
there. Results from this theory were compared, for the cylindrical bending case, with results from 
the exact solution of Pagano 6 for cross-ply laminated plates. The resulting theory, termed a “neo- 
classical” theory, is at least as good as classical theory in. every case and for some cases superior to 
it. Further work was judged to be needed in order to correlate with shear-coupled laminates, also 
treated by Pagano 7 . 


Summary of Work Done During This Period 

Work during this reporting period has continued along two lines: (1) We have continued to 
evaluate the neo-classical plate theory (NCPT) for the shear-coupled laminates and (2) we began 
to explore, with permission from the technical monitor, what kinds of considerations would be 
involved to model plates with nonuniform thickness, which has led to considerable progress toward 
development of higher approximations of our constant thickness plate models. 


Evaluation of NCPT for Shear-Coupled Laminates 

Upon finishing the validation of NCPT for bidirectional plates, we continued to evaluate the 
response of plates with arbitrary stacking sequences. This is a rather challenging problem. Since 
the fiber orientations are not parallel to the axis of each laminate, the influence of shear-coupling 
becomes more evident for this type of laminated plates. When the fiber orientation coincides 
with the principal elastic axes of each laminate no shear-coupling terms exists. For validating the 
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theory, example cases are chosen from Pagano 7 , in which the exact, closed-form three-dimensional 
solutions are available. The results shows that the new theory (NCPT) is more accurate than 
classical laminated plate theory (CPT) when thickness of plate increases. Since NCPT has more 
kinematical variables than CPT, we expect it to be more accurate for thick plates. 

After summarizing the through-the-thickness analysis, we will present the governing equa- 
tions for the global plate analysis with arbitarary stacking sequences. Finally, we will present the 
solutions. 

Three-Dimensional Description Consider a plate of constant thickness h composed of layers, each of 
which is homogeneous and possesses monoclinic material symmetry about its mid-plane; a schematic 
of the plate mid-surface is shown in Fig. 1. Let us introduce Cartesian coordinates Xi so that x a 
denotes lengths along orthogonal straight lines in the mid-surface of the undeformed plate, and 
x 3 = h( is the distance of an arbitrary point to the mid-surface in the undeformed plate, where 
— | < C < Throughout the analysis, Greek indices assume values 1 or 2 ; Latin indices assume 
values 1 , 2 , and 3 ; and repeated indices are summed over their ranges. 

Now, letting bj denote an orthogonal reference triad along the undeformed plate coordinate 
lines, one can express the position vector from a fixed point O to an arbitrary point as 


r{xi,xi ,0 =x a b a + hCb 3 = r(x u x 2 ) + h(h 3 


( 1 ) 


The position vector to the mid-surface is also the average position of points along the normal line, 
at a particular value of xi and x 2 , so that 


r 



( 2 ) 


The angle brackets ( ) are used throughout the paper to denote the integral through the thickness. 
Now, in accordance with Ref. 3 the position vector of any point in the deformed plate is 


R(a;i,£2,C) =R(a:i,a:2) + h(B 3 (x!,x 2 ) 
-\-hWi(xi,x 2 ,C)Bi(xi,x 2 ) 


( 3 ) 


where R is defined as 

R(xi, x 2 ) = (R(zx , x 2 , 0) ^ 

= r(xi,a:2) -I- u(xi,z 2 ) 

and where u is the plate displacement vector, defined as the position vector from a point on the 
undeformed plate mid-surface to the corresponding point on the average surface of the deformed 
plate. The Bj triad is defined so that 


Bi • R,2 = B 2 • R, 1 (5) 

and B3 = Bj x B2 is parallel to (CR)- These definitions give rise to the same kinematical constraints 
on the warping as suggested by Hodges et al . 3 

(Wi) = 0 ((w a ) = 0 (6) 
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We now turn to the strain field, details of which can be found in Ref. 2. First we arrange the 
six strain components into a matrix form so that 


r = [r e 2r, r £ j' 


(7) 


where T e includes the extensional and in-plane shearing strains, and T s and T t contains the trans- 
verse shear and transverse normal strains, respectively. Thus, 


r e = [r„ 2r 12 r 22 j r 
2r* = [2r 13 2r 23 j r 

r £ =r 33 

Next, for notational convenience, we introduce the matrix operators 


( 8 ) 


a i 

dxi 

a 

dX 2 


d e = h 


0 


d t = h 

and generalized strain measures in matrix form given by 


r a_ 

l 

a a 

dXo dxi 


0 afe 


(9) 


£ = 2ei 2 e 22 J r 
K = [Kn 2 ki 2 K 22 J J 
27 = [2713 27 23 J T 


( 10 ) 


where en and e 22 are the plate extensional strain measures, 2ei 2 is the plate in-plane shear strain 
measure, 27 q3 are the plate transverse shear strain measures, K u and K 2 2 are the plate bending 
measures, and 2« 12 = K i2 + K 2 1 is the plate twisting measure. All these measures are functions 
only of Xi and x 2 , and their explicit forms for large deformation are given in Ref. 3. 

Denoting the in-plane warping by 


«>|| = 1^1 V>2 \ T ( 11 ) 

one can now write the Jaumann strain components following Ref. 3 as 


r e = e + ^Kh + d e W\\ 

2 T a = w\ + 27 + d t w 3 (12) 

r t =w ' 3 


where ( )' denotes the derivative with respect to £. 

A similar procedure can be followed for the conjugate stresses so that 

Z e —[Zn Z12 Z22 J T 

Z a = [Z 13 Z 23 \ t (13) 

Zt —Z 33 
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where Z e contains the extensional and in-plane shear stresses while Z 3 and Z t have the transverse 
shear and transverse normal stresses. The stress components may then be written in a matrix form 
as 


Z=[Z e Z, Z t ] T 


(14) 


In light of this, the three-dimensional constitutive law can be expressed as 


Ze 

z a 

Zt 



r D e 

D ea 

D e t 


f r e 

> = 

Dl 

D a 

D a t 


2r, 


Det 

n T 

U 3t 



l r.J 


(15) 


where D e , D es , D et , D a , D st , and D t are 3 x 3, 3 x 2, 3 x 1, 2 x 2, 2 x 1, and lxl matrices, 
respectively. Here, this law is written for directions parallel to plate coordinate axes, which are 
not in general along the material axis. Therefore, the material constants, D's, are the transformed 
values from material axis to the plate axis. 

The plate strain energy per unit area can then be written as 

J = \ (z r r> (16) 


Following Ref. 2, we decompose the strain energy into two positive definite, quadratic forms. Here 
we define the extensional strain energy Jy , and the transverse strain energy J±_ (containing contri- 
butions from both transverse normal and shear strains) as 


Jw = min J 
11 r.,r t 

J± = J - J|| 


(17) 


When the material fiber direction is oriented parallel to the plane of the plate, such that each 
lamina exhibits a monoclinic symmetry, D ea and D at will vanish. In this case the extensional and 
transverse energies can be written in terms of the three-dimensional material properties in the 
following simple form 

2j„ = (rrz7|,r e ) 

2 Jl = (2TJ D a 2Y a + D t (r t + D x Y e f) 

where 

D|| = De - D et D ± D ± = Dt l DZ t 


(18) 

(19) 


This completes the three-dimensional description of the displacement, strain, and stress fields. 
These three-dimensional fields are not really suitable for plate analysis because of the three- 
dimensional warping variables u;*. We now turn to elimination of the warping by dimensional 
reduction. 


Dimensional Reduction In the following sections, we will apply the variational-asymptotical method 
of Ref. 1 for nonhomogeneous, laminated plates in pursuit of a first approximation of the plate strain 
energy. Before doing so, however, it is appropriate to discuss the estimation procedure. First, we 
introduce upper bounds on the in-plane, bending, and transverse shear strain measures e e , £&, and 
e s , respectively, such that 


/ e T c < e e 


^VKTK<e b 


yj 2'y T 2y < e B 


( 20 ) 
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For the first approximation we need only to keep terms in the strain field that are of the order of 
e where 

e e + -f < £ (21) 

This implies that we will have strain energy density of the order fie 2 where fi is of the order of the 
elastic moduli. 

In order to take advantage of the physical aspects of plate deformation, we introduce another 
small parameter h/£, where £ is the smallest constant for which all of the following hold for all 
possible combinations of a and /3 



One may think of £ as the wavelength of the deformation pattern. 

Rather than write out complete expressions for the strain field, we will only write the terms 
needed for the first approximation. By consideration of the above set of estimation parameters, the 
strain expression can be approximated as 


T c = e + QKh 

2r„ = 27 + tu(| (23) 

r t = wa 

Thus, the warping is not present in extensional energy, and the only part of the strain energy 
remaining to be minimized is the transverse energy 

2Ji=(AK + Z?x (e + C Kh)] 2 ) 

/ \ (24) 

+ ((27 + m[|) T A(27 + w[|)) 


The variational-asymptotical method calls for the minimization of this strain energy expression 
with respect to the warping, with the constraints given in Eq. (6), resulting in 


w 3 =D_ u€ + DsjzKh 


w ll = 


~(Hs 1 - 8f/, 2 )«CAi) - 8(C H s2 ))~ l - / 2 C 


27 


(25) 


in which / 2 is the 2x2 identity matrix and 

D' Ll =-D ± D'j_ 2 = -(D _ l 

C 2 (26) 

K 1 = D,- 1 H' 32 = yA -1 

These expressions for the warping are determined uniquely by imposing the continuity of D± a and 
H sa between the layers together with constraints Eq. (6) so that 

< Al«> = ( H . a ) = 0 (27) 
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for out-of-plane and in-plane warping, respectively. These constraints guarantee that warping 
functions are continuous. However, strain and stress, which are functions of the derivatives of the 
warping, may or may not be continuous. 

The strain energy per unit area of the plate is obtained by substituting the warping from 
Eq. (25) into the strain energy, yielding 

J = i [e T (Ae + 2 BK) + K T DK + 2 7 T G2 7 ] 

where 

A = h(D {] ) B = h 2 (CZ)||) 

D = h 3 (C 2 Z?|,) G = h(g T D a g) 

and where 

9 = ~ (i - 4C 2 ) ».) - S(C^))- 1 

Note that D\\ corresponds to Q, the well known transformed reduced stiffness matrix from classical 
laminated plate theory 8 . It is possible to define the force, moment and transverse shear stress 
resultants N, M, and Q respectively, as 


(28) 

(29) 

(30) 


rp 

N = h(Z e ) = 


Q = h(Z.) = 



(31) 


Then, based on the strain energy, the plate constitutive law can be expressed as 


m 

' A 

B 

O' 


< M > = 

B t 

D 

0 

\ * 

UJ 

0 

0 

G 

UtJ 


Note that transverse normal stress is zero in this theory, and thus, we should not expect the 
transverse normal strain to be very accurate from Eqs. (12) unless we extend the theory to higher 
approximations. Similarly, we do not expect the transverse shear stress and strain to be very 
accurate with the present theory, since some important contributions to its detailed variation are 
associated with h/i corrections to the energy 1 . These quantities can be obtained from integrating 
the three-dimensional equilibrium equations through the thickness to get the transverse shear and 
normal stresses and applying the three-dimensional constitutive law to get the transverse shear and 
normal strains. 

This concludes the dimensional reduction. The global deformation equations 3 , along with 
these plate constitutive equations, comprise what we term the neo-classical theory. (Note that the 
displacement shift mentioned in Refs. 2 and 3 is not necessary with this theory, and thus 2 7 * = 2 7 .) 
With the warping known in terms of e, K, and 2 7 , which in turn are known through solution of 
the global deformation problem, it is now possible to evaluate three-dimensional approximations of 
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displacement, strain, and stress fields. For the purpose of validating the stiffness model and field 
relations, however, only a specialized version of the global deformation analysis is undertaken here. 


Linear Plate Equations The global deformation equations that correspond to the above strain 
energy function were developed in Ref. 3. Since e, K , and 27 are nonlinear functions of the 
displacement and rotation variables, this theory is applicable to large deformation of plates, and 
these can now be obtained by solving a specific problem. Here, rather than repeat the entire 
formulation, we will specialize it for linear, cylindrical bending problems, which we will use below 
for validation of the dimensional reduction scheme of the previous section. 

Kinematical equations from Ref. 3 for the linear case are given as 


til = « 1,1 
2 ei 2 = « 1,2 + ^ 2,1 
£22 = « 2,2 
Kll = # 1,1 

2 «i 2 = # 1,2 + $ 2,1 

K22 = $ 2,2 

2713 = $1 + « 3,1 
2723 = $2 + u 3,2 

where n, = u • b t - and 0 a = B 3 • b Q . 

Similarly, equilibrium equations are 

Nn.i + Ni2,2 + h — 0 
1^12,1 + 1^22,2 + /z = 0 

Ql,l + ^ 2,2 + h = 0 

A^n.i + Mi2,2 — Qi = 0 
■Mi2,l + M22,2 — Q 2 = 0 

where fi are the applied loads per unit area of the plate. 


(33) 


(34) 


Cylindrical Bending Analysis Consider a plate of length L along x\ and infinite width in the a ; 2 
direction shown in Fig. 2. All derivatives with respect to 22 are zero which causes many of the 
variables in the above equilibrium, constitutive, and kinematical equations to drop out. 

In the case of cylindrical bending, the plate is subject to sinusoidal surface loading of the form 

/a(s 1 ) =Po sin(pii) (35) 


where V — z anc ® fa = 0- For the three-dimensional analysis, we assume the loading to be 
imposed in the form of an upper surface traction. 

We consider a simply supported plate, the boundary conditions of which are 

u 3 (0) = u 3 (L) = 0 ««,i(0)=u a ,,(L) = 0 $«,i(0)=$ a ,i(L) = 0 (36) 
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The governing system of equations can be put into a matrix form by defining generalized coordi- 
nates, strain measures, stress and moment results, and loading. Let 


<?= [til U2 0\ &2 Us\ T 

(37) 

£=[cn 2ei2 «u 2 k 12 2713 2723 J T 

(38) 

F= [N n N 12 Mu M 12 Qi £?2j T 

(39) 

h 

I 

0 

0 

0 

0 

II 

(40) 


Eqs. (33) can be written as 


e = V q (41) 

where 

9 0 0 0 0 

0 3 0 0 0 

0 0 3 0 0 

0 0 0 3 0 

0 0 10 3 
0 0 0 1 0 

in which 3 denotes the derivative with respect to Xi- Eqs. (34) can be written as 

f = £ F (43) 

where 

'3 0 0 0 0 O' 

0 3 0 0 0 0 

£ = 0 0 3 0 -1 0 = 3* t - k T (44) 

0 0 0 3 0 -1 

.0 0 1 0 3 0 . 

Now plate constitutive equation can be written as 

F = K*e (45) 

where 

A* B* O' 

K* = B” D* 0 (46) 

0 0 G 

Here the starred matrices are 2x2 sub-matrices consisting of the first 2 rows and columns of the 
corresponding matrices in Eqs. (32), resulting in K* being a 6 x 6 matrix. So the given problem 



rsoooo- 
0 3 0 0 0 

0 0 3 0 0 

" 0 0 0 3 0 

0 0 0 0 3 

LO 0 0 0 0. 


'0 0 0 0 0 ' 

0 0 0 0 0 

0 0 0 0 0 

+ 0 0 0 0 0 
0 0 10 0 
.0 0 0 1 0 . 


= 3* + k (42) 
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yields 17 equations in terms of 17 unknowns. Combining Eqs. (41), (43), and (45), the equilibrium 
equations can be written as 

(d* r - K * ( d * +k)=f (47) 


which can be put into the form 


where 


M n = 


M c q , a + C c <7 iX - K c q = f 


•^11 

7^12 

1 

Bl 2 

0 - 




^12 

■^22 

B{ 2 

b 22 

0 


A* 

B* 

0 ‘ 

B'n 

Bit 

DU 

Dh 

0 

= 

B* 

D* 

0 

Bh 

Bh 

D\ 2 

Dh 

0 


0 

0 

G„ 

. 0 

0 

0 

0 

Gn- 





K, = 




■0 0 

0 

0 


0 - 





0 0 

0 

0 


0 



nr 

'c — 


0 0 

0 

0 

- 

Gn 





0 0 

0 

0 

- 

G 12 





.0 0 

Gn 

G 12 

0 . 



0 

0 

0 

0 

01 





0 

0 0 

0 

0 


'0 

0 

o' 

0 

0 Gn 

G 12 

0 

- 

0 ( 

jr 

0 

0 

0 

G 12 

G 22 

0 


0 0 

0. 

.0 

0 

0 

0 

0. 






(48) 


(49) 


(50) 


(51) 


The above equation can be solved directly with associated boundary conditions. But let us take 
advantage of the sparcity of matrices in the equation. We cam break down this equation as follows 


A* tt,n +B* 9,u = 0 


(52) 


B* ti,ii + D* 6, ii — [Gn G 12 J T U3,i — G 9 = 0 

Gn U3,n + [Gh G 12J 6 = — /3 


(53) 

(54) 


when u and 0 are arranged as 


u=[ui u 2 \ T 9 = [61 0 2 J 7 (55) 

Since / 3 is a simple sinusoidal function given in Eq. (35), 9 can be uncoupled by putting Eqs. (52) 
and (54) into Eq. (53) yielding 


Ml 9,ii + K 9 = SI j cos (pan) 


( 56 ) 


where 


AC = D* - B* A* ~ X B' 


K = 



fc = L 1 
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Since M * is symmetric, positive definite and K* symmetric, Eq. (56) can be decoupled as is. Note 
that when plate is symmetrically layered up B* matrix vanishes, resulting M* = D*. 

Now, introducing the coordinate transformation 

0 = 4> y ( 57 ) 


where 4> is a eigenvector matrix, we diagonalize Eq. (56) which gives 

M d y,u+K d y = cos (px x ) (58) 

where 

M d = <t> T M:<j> K d = <f) T K*<f) f d = <t> T f' c (59) 

By solving above decoupled equations for y, we can recover original solution 6 using Eq. (57). It 
is easy to shown that all integration constants produced in the calculation vanish for cylindrical 
bending. In other words, we need only to obtain particular solution for Eq. (56) or Eq. (58). Note 
that since effects of boundary layer belong to higher order for this type of first approximation, we 
can neglect a homogeneneous solution in the first place. 

After some calculation, 6 can be found. Also u 3 and it are determined by putting 0 into 
Eqs. (52) and Eqs. (54). Rotations, in-plane displacements, and out-of-palne displacement are 
obtained as follows 

0 = 4> [K d - p 2 M d ] 1 f d ^ cos (pxi) (60) 

it = -p 2 M d ]~ l f d cos (px 1 ) (61) 

it 3 = ~~ (l - LGn G 12 J <(> [K d -p 2 M d ]~ l psinfazi) (62) 

To recover three-dimensional fields, one can substitute values obtained for e, K , and 27 by 
taking derivatives of above expressions into Eqs. (25), (23), and (15), to get warping, strain, and 
stress, respectively. Note that all the strain measures can be calculated to order e by using Eq. (23); 
however, r e can be determined to order he/i by use of Eq. ( 12 ). To avoid inconsistency, one should 
use Eq. (23) for all strain components when calculating the stress components. However, NCPT 
is more accurate than CPT since it picked up the general strain measures which do not appear in 
CPT. In other words, the extra strain measures reflect the effects of transverse shear deformatiom, 
resulting in a better solution: 

For the purpose of comparing the displacement field with that of linear, three-dimensional 
elasticity, we need to determine the bi measure numbers of the displacement field 

Zi = bj • (63) 

With the warping known from Eqs. (25) in terms of e, K, and 27 , and the triad Bj given as 

Hi = bi -1 — ^^2 — 0i b 3 

B 2 = - ^-bj +b 2 - 0 2 b 3 (64) 

B 3 = $ibi — 02^>2 T b 3 
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The functions Zi can now be expressed in terms of global deformation variables. The result is that 

zi = ui+ h^di + hw i 

z 2 = U2 - h(02 + hw 2 (65) 

z 3 = u 3 + hw 3 

where Ui, & a are obtainable in Eqs. (60), (61), and (62), and Wi in Eq. (25). 


Results The intent of this section is to compare results from the exact solution with results from 
the present theory. The exact three-dimensional elasticity solution for cylindrical bending of angle 
ply laminated plates was obtained by Pagano 7 , and the results presented herein labelled as “exact” 
generated from his equations. The material properties are 7 


El = 25 x 10 6 psi Er = 10 6 psi 

Glt = 0.5 x 10 6 psi Gtt = 0.2 x 10 6 psi (66) 

vlt — vtt = 0.25 

where L signifies the direction parallel to the fibers and T the transverse direction. We evaluated 
all quantities obtainable from the neo-classical theory for this problem and compared those results 
with the exact solution for various values of L/h. We have chosen to present most of our results 
for L/h values of 10 and 4, which represent relativatively thin and thick plates, respectively. 


For plotting the displacement, strain, and stress distributions the following normalized param- 
eters are used 

E T h 2 z a (x u O ^ 100ETh 3 z 3 (x u O 

Za ~ L 3 p 0 23 L*p 0 

— h 2 Z a p{x 1,C) -7y hZ a3 (x\,£) fan\ 

Za * = L 2 po Za3 = " Lpo (6?) 

_ E L h? r^.fo.c) 
u L 2 p 0 


The through- the-thickness distributions of in-plane displacement, transverse shear stress and 
strain are evaluated at xi = 0; while the distributions of out-of-plane displacement, in-plane stress 
and strain, transverse normal stress and strain are evaluated at X\ = L/2. In all the results below, 
solid lines represent the exact solution, while dashed lines represent the present neo-classical plate 
theory (NCPT) results. Results from classical laminated plate theory (CPT), when distinct from 
NCPT results, are shown with long and short dashes. 

The distributions of normalized displacements, strains, and stresses are obtained for the fol- 
lowing three plates, in which each layer has the same thickness: 

[15°J 

[15°/ ~ 15°) 

[307-30°] aym 

Results for these plates are presented in Figs. 3 - 19. Figs. 3-5 show the transverse displacement 
for these three configurations, respectively. Note the very close agreement between NCPT and the 
exact solution. NCPT is just slightly more flexible than the exact solution; it provides significant 
improvement over CPT as L/h becomes small, until L/h of approximately 10. 
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Figs. 6 and 7 show normalized out-of-plane displacement (z 3 ), in-plane displacement (zi), in- 
plane strain (rn) and shear strain (2Fi2) for (15°j plate from top to bottom in the figures; Figs. 10 
and 11 show the same quantities for (15°/ — 15°] plate. Figs. 8 and 9 show normalized in-plane 
stress (Zn), transverse shear stress ( Zi 3 ) and strain (2Ti 3 ), transverse normal strain (r 33 ) for [15°] 
plate; Figs. 12 and 13 show the same quantities for [15°/ — 15°] plate. 

Figs. 14 - 19 show the results for [30°/ — 30°]^^ plate. Figs. 14 and 15 show normalized 
out-of-plane displacement (z 3 ), in-plane displacements (zi, z 2 ), in-plane strain (r n ). Figs. 16 and 
17 show normalized in-plane shear strain (2r 12 ), in-plane stress (Zn) and shear stress (Z 12 ), and 
transverse shear stress ( Z 13 ). Figs. 18 and 19 show normalized transverse shear stress (Z 23 ) and 
strain (2r 23 ), transverse normal stress (Z 33 ) and strain (r 33 ). 

Transverse in-plane stresses (Z22) are not shown since these have the same pattern as the cor- 
responding in-plane stresses (Zn). Note that transverse shear stresses in the figures are obtained 
from integrating the equilibrium equations of elasticity. Also note that transverse in-plane dis- 
placement (z 2 ), in-plane shear strain and stress (2 Fj 2 , Z 12 ) and transverse shear strain and stress 
(2r 23 , Z23) will differ from those of Pagano 7 . This stems from the different sign convention used 
for plotting in this work, resulting in the same magnitude but with opposite sign. Some of these 
quantities are also reported in graphical form in Ref. 9. 

The in-plane displacements and strains from NCPT shown in Figs. 6, 7, 10, 11, 14 and 15 are 
much closer to the exact solution than those from CPT. One reason the in-plane strain and shear 
strain are so accurate is that F e is calculated to 0(he/t) using Eq. (12). The in-plane stress and 
shear stress are essentially identical to that of CPT shown in Figs. 8, 9, 12, and 13 for the 1- and 
2-layer cases. On the other hand, NCPT is better than CPT shown in Figs. 16 and 17 for the 
4-layer case. This is because that NCPT has more kinematical varivables than CPT. Even though 
the difference is invisible for 1- and 2-layer cases, twisting curvature, which is a one of generalized 
strain measures used in NCPT, plays an important role for this specific example. 

It is interesting to note that in Fig. 19, transverse shear stress (Z 23 ) from NCPT has sign 
changes around at the middle surface, while results from the exact solution are positive. Based 
on results obtained to date, we believe that an extension of the theory of Ref. 3 to a higher 
approximation still needs to be developed in order to improve the correlation of this and other 
three-dimensional field variables with the exact solution. However, interaction terms are present 
in the energy when a higher approximation is attempted, similar to those of Ref. 2; these must 
be killed, but the means for doing so were unknown until recently. This means was discovered in 
conjunction with the analysis of nonuniform thickness plates. 


Higher Approximations and Plates with Nonuniform Thickness 

During the last two months, we have been attempting to extend our methodology to deal with 
plates, the thickness of which varies spatially over the plate domain so that h = h(x i,x 2 ). In the 
process of doing this for isotropic plates, an additional degree of freedom for the normal line element 
was introduced on physical grounds, within the context of VAM. This degree of freedom involves 
the contraction of the normal line element. Although the contraction was not zero in the original 
theories, it was a “reactive” quantity. That is, it does not appear in the strain energy function 
explicitly, and it can be calculated in terms of other degrees of freedom of the normal line element. 
When the thickness is allowed to vary spatially over the plate, this degree of freedom appears 
explicitly in the reduced, two-dimensional strain energy function and must be regarded as an 
independent quantity. This new degree of freedom improved the correlation with three-dimensional 
results relative to classical theory, but the accuracy was still judged inadequate. Indeed, it was not 
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possible to show that the additional warping was of a higher order than this degree of freedom. 
Moreover, the motivation for introducing this degree of freedom was not at all systematic. 

Dr. Sutyrin then discovered a powerful principle for identifying valid degrees of freedom in the 
dimensional reduction process. (This principle applies for development of beam, plate, and shell 
theories.) There exists an eigenvalue problem, the eigenvectors of which are the valid degrees of 
freedom for the structural member. They are valid in the sense that (1) they satisfy the appropriate 
interfacial conditions; (2) they guarantee that the above-mentioned interaction terms vanish; (3) 
they guarantee that any additional warping is of a higher order. They are not simple polynomials, 
in general. The degrees of freedom in classical theory, it turns out, are the eigenvectors associated 
with the zero eigenvalues for this eigenproblem. Sutyrin’s eigen-principle should allow us to extend 
the earlier laminated plate theories to higher approximations in a rational manner. 

If plates of constant thickness must be joined to plates of varying thickness, then both models 
must contain the same degrees of freedom. This means that the constant thickness model must 
be taken to a higher-order approximation so as to bring the same, or at least analogous degrees of 
freedom into the strain energies of both plates. The mathematics associated with this operation are 
quite difficult, involving functional analysis as well as the construction of non- trivial transformations 
to rid the strain energy of unneeded degrees of freedom. Unanswered questions also remain con- 
cerning the role of edge-zone behavior and penetration of the influence of self-equilibrated tractions 
on the plate edges. 

To illustrate this, let us write the Jaumann strain components as 

T e = e + £Kh + d e Wft 

2T a = io[j + d t w 3 (68) 

T t = w' 3 

The resulting warping will be different from the result of Eqs. (25), because 2y has been set to 
equal to zero and its effect absorbed into w. It turns out that the weirping variable can be written 
as 

Wi(x 1,X 2 ,0 = J2^(x 1> x 2 )ip ik (0 +Vi(zi,X2,0 (69) 

fc 

where v, represents additional warping and q, the degrees of freedom associated with the ip functions. 
The ip’s can be identified by means of the VAM procedure and Sutyrin’s eigen-principle. The 
resulting two-dimensional strain energy function will contain the q functions and their derivatives. 
The process of identifying these “degrees of freedom” for the plate, as presented in the literature, 
is not straightforward. We have made some important progress recently in this identification 
procedure. 

We consider a homogeneous, isotropic plate with thickness depending on only one in-plane 
coordinate and undergoing a stretching load; see Fig. 20. In order to solve this problem we have 
to introduce new degrees of freedom. It is clear that the average bending of the plate and rotating 
of normal line element are equal to zero here. Nevertheless there is a non-trivial deformation of 
normal line elements, which can be observed in the elasticity solution (which we approximated by 
a finite element solution). This deformation is shown schematically in Fig. 21. There is no known 
systematic way to introduce such deformations into classical theory or into Reissner’s theory, which 
takes into account only the transverse shear deformation. 

In order to introduce some new degrees of freedom in our analysis, we take the warping to be 
of the form of Eq. (69) where ip is the shape function of a degree of freedom and q is the new degree 
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of freedom itself. The VAM procedure points to the best choice of the shape function ip as an 
eigenvector of a natural eigenvalue problem for normal line element (a one-dimensional problem). 

For example, the two degrees of freedom corresponding to the smallest nonzero eigenvalues 
can be represented as 


w a = 9a(xi , x 2 ) cos(2ttC) + V a 

(70) 

w 3 = e(xi,x 2 )sm(TrO + v 3 

(71) 


for the homogeneous isotropic case. The degree of freedom g x ( g for short) is responsible for in- 
plane deformation such as shown in Fig. 21, while the degree of freedom e describes contraction of 
the normal line element. The variable v models any additional warping, which is of a higher order 
and may be ignored in this example. 

Fig. 22 shows comparison between g and e obtained from a 3-D solution (dashed line) and g 
and e from our new theory (solid line) . 

This approach has a natural generalization for dynamics problems. The above eigenvalues 
become eigenfrequencies in that context. Also, it is not difficult to obtain a suitable dynamical 
equation for each degree of freedom. In general, one would need to apply the theory to a wide class 
of problems and loading conditions in order to understand the appropriate number of degrees of 
freedom to be retained. 

We expect to develop the dimensional reduction with the extensive use of computerized sym- 
bolic manipulation and the expertise of Dr. Sutyrin. As a fallback position, should a symbolic 
dimensional reduction not be possible, we can solve the minimization problems via one-dimensional 
finite elements (through-the- thickness). Both means are quite efficient for laminated plate prob- 
lems. The through-the-thickness analysis is only done once for a given lay-up. The results for 
the elastic constants are used as input to the plate (two-dimensional) analysis, and the influence 
functions are used to recover approximations for the three-dimensional field variables once the 
two-dimensional problem is solved. 


Future Work 

In the balance of the grant we intend to apply Sutyrin’s eigen-principle to develop a refined 
theory for laminated plates of constant thickness. If there is time we will make this theory capa- 
ble of modeling the most general type of nonhomogeneous, anisotropic plate, subject only to the 
restrictions that the strain is small and that neither geometry nor properties vary with X\ and x 2 - 
Development of an interior global deformation analysis for constant thickness anisotropic plates 
would then be possible, but it appears that this should be addressed only after the degrees of free- 
dom for a general analysis are identified. This means that the boundary energy and nonuniform 
thickness problems should be addressed first. We are seeking to obtain funding from NASA to 
carry this out. 
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Fig. 12: Distributions of Normalized Quantities for 
[15°/ — 15°] Plate when L/h = 10 
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Fig. 13: Distributions of Normalized Quantities for 
[15°/ - 15°] Plate when L/h = 4 








Fig. 16: Distributions of Normalized Quantities for Fig. 17: Distributions of Normalized Quantities for 
[30°/ - 30°1 Plate when L/h = 10 [30°/ - 30°j Plate when L/h = 4 
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Fig. 20: Tapered plate subjected to large in-plane force 



Fig. 21: Exaggerated deformation of normal line elements for tapered plate 
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